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The density of states is calculated for a SU(2) and a compact U(l) lattice gauge theory using a 
modified version of the Wang-Landau algorithm. We find that the density of states of the SU(2) 
gauge theory can be reliably calculated over a range of 120,000 orders of magnitude for lattice sizes 
as big as 20 4 . We demonstrate the potential of the algorithm by reproducing the SU(2) average 
action, its specific heat and the critical couplings of the weak first order transition in U(l). 
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Monte-Carlo simulations [T] of the theory discretised 
on a Euclidean space-time lattice [2] currently provide 
the most successful approach to calculations from first 
principles in asymptotically free gauge theories in the 
energy domain in which the coupling is of order one. 
Although this strategy is successful for computations of 
observables that can be expressed as a vacuum expecta- 
tion value (vev) on a theory with a semi-positive definite 
path integral measure, when the observable is not a vev 
(e.g., the free energy, which is related to the logarithm 
of a partition function) or the path-integral measure is 
not semi-positive (like in QCD at finite density), Monte- 
Carlo algorithms are either unsuitable or very inefficient. 

An alternative numerical approach to Lattice Gauge 
Theories potentially free from those limitations is based 
on the density of states. Let us consider a quantum field 
theory with action pS[<j)\, with (3 the inverse coupling. 
For this theory, the path integral in Euclidean space-time 
is given by 



Z 



V(j>{x) 



(i) 



where (T>cj)(x)) means that the integral has to be per- 
formed over all allowed configurations of the field </>. 
Defining the density of states p(E) as 

p{E) = J V<t>(x) 5(S[4>] - E) , (2) 

the path integral can be rewritten as 

Z I p(E)e 0E dE, (3) 



and the vev of an observable O(E) becomes 

(0) = i J p{E)0{E)^ E dE . (4) 

If the density of states is known, the path integral and 
(O) can be obtained by computing numerically or ana- 
lytically respectively the integral ^ and Q. 

An efficient algorithm for computing p(E) in systems 
with discrete energy levels has been proposed by Wang 



and Landau in To date, the method has found vari- 
ous applications in Statistical Mechanics, some of which 
have produced remarkable results that can not be ob- 
tained with a direct Monte-Carlo approach (see e.g. [3] 
for a recent example). Despite its popularity in Statis- 
tical Mechanics, the Wang-Landau algorithm has found 
only limited applications in Lattice Gauge Theory [5J [B] . 
In fact, the sampling of a continuous density of states 
with a straightforward generalisation of the method given 
in [3] turns out to be problematic [8] . In this work, 
we propose a new method for determining a continuous 
density of states and we apply it to calculate the density 
of states in SU(2) and U(l) on the lattice. 

Throughout this paper we adopt the lattice regular- 
isation, which leaves us with a A 4 cubic lattice as the 
discretisation of the Euclidean space-time. The dynam- 
ical degrees of freedom of the SU(A C ) gauge theory are 
represented by the matrices U^x) £ SU(A C ), which are 
associated with the links of the lattice. We are using the 
so-called Wilson action, i.e., 
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E 

fJ 1 >l',X 



1 



Retr 



U^x)U v {x + p)Ul{x + v)Ui{x) 
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stressing however that our approach is not limited to this 
particular action, but can handle e.g. improved actions 
equally well. 

In order to present our novel type of numerical algo- 
rithm to calculate the density of states, we will assume 
that In p{E) is well approximated by piecewise linear 
functions. It will indeed turn out below that In p(E) 
is a remarkable smooth function of E. 

Let us consider the energy interval [E ,E + 5E] for 
which we approximately write 



p(E) = p(E ) cxp{a(£ ) (£-£„)} 



(6) 



for Eq < E < Eq + SE. Our goal will be to calculate the 
coefficients a(E ), which can be considered as derivatives 
of the density of states: 



a(E ) 
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FIG. 1. The thermalisation history for a SU(2) gauge theory 
for lattice sizes 10 4 . . . 20 4 . 



FIG. 2. The statistical error for the estimate of a(Eo) for 
lattice sizes 10 4 ...20 4 . 



The strategy to obtain these coefficients is based upon the 
truncated and re-weighted expectation values defined by 

«/(£)» (a) = ^JdE f(E) p(E)6 [EotSE] e- aE , (8) 

N = J dEp(E)6 [Eo>SE] e~ aE , (9) 

f 1 for E < E < E + 5E, f . 
0[B ,SE] - | elsewhere . 

If the energy interval is small enough, i.e., if © is a 
good approximation, we should be able to choose a to 
compensate a(Eo). This would leave us with a flat energy 
histogram and with 

((E))(a) = E Q + — , for a = a(E ). (11) 

Assume now that a n is an approximation for a(E Q ) such 
that x = [a(E ) - oj SE < 1. Defining AE := E — E — 
SE/2, we then find using ^ 

SE 2 

((AE))(a n ) = — [a(E ) - a n ] + 0(x 3 SE) . (12) 

Ignoring the higher order correction and solving for 
a(E ), we obtain a better approximation a n+1 : 
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(13) 



The central idea is to iterate the latter equation until 

({AE))( aoo ) = =>• aoo =a(Eb) , 



where we have used (11). We point out that the trun- 



means of Monte-Carlo methods. To this aim, we insert 
^ into ([8} to obtain: 

((f(E)))(a) = i / VU, f(s[U}) e- aS M ,(14) 
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The subscript of the integral indicates that updates of 
configurations the action of which falls outside the de- 
sired energy interval are discarded. There are many 
Monte-Carlo techniques to estimate the truncated expec- 
tation value in ( [14] ), the Metropolis algorithm and the 
Heat-Bath approach being the two most obvious choices. 
We have tested both techniques and found that our 
method for estimating o,(Eq) is robust. The numerical 
results shown below have been obtained by an adapted 
Heat-Bath algorithm with a 100% acceptance rate (de- 
tails of the algorithm will be published in a forthcoming 
paper). 

Let us now consider the SU(2) gauge theory to il- 
lustrate our approach in practice. If iV 4 is the num- 
ber of lattice points, the maximal action is given by 
B max = 6iV 4 . We here consider the energy interval 
I := [E ,E + SE] = [0.650, 0.651] 6N 4 . The first task 
is to generate a lattice configuration {U^} the action of 
which falls into the energy interval /. For this purpose, 
we start with a "cold" configuration Uu(x) = 1, and up- 
date the configuration forcing it to reach the desired en- 
ergy interval. We then pick a start value for the iteration 
(13), which has been oq = —2 in this preliminary study. 



cated expectation values can be easily estimated by 



We perform 25 energy restricted Monte-Carlo sweeps at 
do (see (14)), where each sweep consists of iV 4 updates 
of randomly chosen individual links. 

In order to evaluate the next Oi the expectation value 
((AE)) is evaluated using the energy restricted Monte- 
Carlo method (see (14)). For this, we have used 384 
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FIG. 3. The estimates for a(J3 ) for S = 0.650 x 6iV 4 as a 
function of the lattice size. 
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FIG. 4. The logarithm (base 10) of the density of states for 
the SU(2) gauge theory using a 10 4 lattice. 



measurements divided in 48 independent runs each con- 
tributing 8 Monte-Carlo sweeps (these calculation are 
performed on the HPC computing facilities at the Ply- 
mouth University). The corresponding estimator is then 
used to obtain an improved value a%. This procedure 
is reiterated n times, n > 1, until the value of a starts 
to fluctuate around a central value. The thermalisation 
history is shown in figure [T] for small lattice sizes such 
as 10 4 , a thermalised state is reached after 10 iterations 
while for our biggest lattice 20 4 roughly 80 iterations are 
necessary to reach an equilibrium. To keep control of the 
autocorrelation in the determination of the solution of 
the iterative procedure we have evaluated the integrated 
autocorrelation time (T int ) of ((AE)). In particular the 
measure of Ti nt for the highest energy gap yields a value 



FIG. 5. Average plaquette for a SU(2) gauge theory on a 10 4 
lattice obtained by means of the density of states and local- 
hybrid Monte-Carlo. Also shown is the specific heat x(P)i 
Eq. p|. 



always smaller that two steps for each of our volumes. 

Having control of the autocorrelation time allows us to 
reliably define a statistical error of ((AE)) which directly 
feeds into the uncertainty for a n+ \ (see (13l). Rather 



than to spend all numerical resources to obtain a high- 
precision estimate for ((AE)) we found is advantageous 
to feed the more noisy estimator into the iteration ( 13 ) 



and to average the a n values of the resulting sequence. 
The standard error of a n for an average over a bin of 10 
iterations after thermalisation is shown in figure [2] We 
roughly find that the error decreases like 1 / VV where V 
is the lattice volume. The lack of autocorrelation reflects 
in the good scaling of the error with the volume show- 
ing the efficiency of the algorithm also for large volumes. 
In particular, this observation is true even when study- 
ing energy intervals for which we would normally expect 
strong effects in autocorrelation due to critical slowing 
down (for example for 0.850 < E/Exaax < 0.851 and 
V = 20 4 we find r mt = 1.8(1)). 

For the determination of a(Eo), 187 iterations have 
been performed for thermalisation and 312 further iter- 
ations were carried out to estimate cl(Eq). Our findings 
as a function of the lattice size are shown in figure [3j 



Once cl(Eq) has been obtained for all energies Eq 



i x 



SE (we here only consider positive energies), the density 
of states p(E) can be easily constructed from ([6]): 



p(e)=h< 



,a{E' )SE 



i = l 



exp{a(£ ) (16) 



for Eq < E < E% +1 . Thereby, we have normalised the 
density of states such that p(E = 0) = 1. Our numerical 
result is shown in figure |4j In order to estimate any in- 
fluence of the discretisation error, we have calculated the 
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FIG. 6. The probability density Pp{E) for a compact U(l) 
gauge theory at critical coupling for several lattice sizes iV 4 . 



density of states by splitting the energy interval [0, E max ] 
into 1000 and 5000 energy intervals. Both curves fall on 
top of each other in figure |4| As a proof of concept that 
our numerical approach does yield high precision expec- 
tation values, we have calculated the average plaquette 
{E)/E max using Q. As expected, only a small energy 
window with a(E) « j3 significantly contributes to the 
expectation value. Care has been taken to handle poten- 
tially large numbers. We have compared our result with 
that from a standard method using local-hybrid Monte- 
Carlo. A very good agreement is observed. An observable 
which is generically difficult to estimate due to cancella- 
tions is the specific heat, which we define by 

X{P) = g^i (V) - ^) 2 )' ( 17 ) 

where the expectation values are obtained by means of 
Q. Our numerical findings for x are also shown in fig- 
ure [5] We have checked for a few (3 values that our result 
agrees with that obtained by standard methods. 

We have finally tested our approach for the compact 
U(l) gauge with Wilson action (p5]). Here, the links are 
U(l) group elements, i.e., U^(x) — exp{z6* M (x)} with 
@n( x ) = — 7T . . . 7r being the dynamical degrees of free- 



dom featuring in the functional integral with a constant 
measure. By means of a large scale investigation on 
the basis of the Borgs-Kotecky finite size scaling anal- 
ysis, it has been finally established in [5] that compact 
U(l) possesses a weak first-order phase transition at 
P = /3 C rj 1.0111331(21) (in the infinite volume limit). 
An unmistakable sign for a first order transition is the 
characteristic double-peak structure in the action proba- 
bility density, i.e., 

Pp{E) = p(E) exp{/3£} , (18) 

for p — > p c . It turns out that this double-peak structure 
is very sensitive to variations of p allowing a high preci- 
sion determination of P c at finite volume, i.e., the critical 
coupling for which the peaks are of equal height. Note 
that we have normalised P/s{E) such that its maximum 
value equals one. The critical couplings /3 C , listed in the 
graph, are in good agreement with those from the large 
scale study [9]. 

In conclusions, we have developed a modified version 
of the Wang-Landau algorithm suitable for theories with 
continuous degrees of freedom. We have shown that the 
density of states for a SU(2) gauge theory can be calcu- 
lated over a range of 120,000 orders of magnitudes even 
for a lattice as large as 20 4 . Our approach reliably re- 
produces the critical couplings of the weak first order 
transition of the compact U(l) gauge theory. Using the 
Cabibbo Marinari method [TU], our approach can be gen- 
eralised to SU(iV c ) Yang-Mills theories. Quantities of 
interests which are earmarked for our approach are ther- 
modynamic potentials |11) . vortex free energies |12j and 
electric fluxes for the study of the mass-gap and confine- 
ment |13j . Finally, we point out that the statistical error 
for expectation values obtained by the density of states 
method can be obtained by the bootstrap technique. A 
careful investigation of the statistical and possible sys- 
tematic errors (from which our results seem to be free) 
will be reported elsewhere. 
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